Erythroid clean

# Seurat 
library(Seurat)

# Data 
library(dplyr)

# Plotting
library(ggplot2)
library(RColorBrewer)
library(patchwork)
library(gridExtra)
library(grid)
library(viridis)
library(ComplexHeatmap)

Attaching SeuratObject


Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union



Attaching package: ‘gridExtra’


The following object is masked from ‘package:dplyr’:

    combine


Loading required package: viridisLite

========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================

# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")

# Source files
source("plotting_global.R")

Parameter settings

# Seurat files 
so_file <- "data/object/seurat.rds"
so_out_file <- "data/object/seurat_clean.rds"

# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

Import Seurat object

so <- readRDS(so_file)

# Create clean Seurat object
so <- CreateSeuratObject(counts = GetAssayData(so, assay = "RNA", slot = "counts"), meta.data = dplyr::select(so@meta.data, -rna_snn_res.0.8))

Subset by sample group

so_l <- SplitObject(so, split.by = "treatment")

so_l <- lapply(so_l, NormalizeData)

so_l <- lapply(so_l, FindVariableFeatures, nfeatures = 2000)

cc_kowalczyk <- read.csv("cc_kowalczyk.csv")
cc_kowalczyk <- cc_kowalczyk[cc_kowalczyk$sig_population >= 5, ]

variable_features_reduce <- function(so, genes) {
    variable_features <- VariableFeatures(so)
    VariableFeatures(so) <- variable_features[!variable_features %in% genes]
    return(so)
}
so_l <- lapply(so_l, variable_features_reduce, genes = cc_kowalczyk$gene)

so_l <- lapply(so_l, ScaleData, vars.to.regress = c("nCount_RNA"))

Regressing out nCount_RNA

Centering and scaling data matrix

Regressing out nCount_RNA

Centering and scaling data matrix

options(warn = -1)
so_l <- lapply(so_l, RunPCA, npcs = 100, verbose = FALSE)
so_l <- lapply(so_l, FindNeighbors, dims = 1:30, verbose = FALSE)
so_l <- lapply(so_l, FindClusters, verbose = FALSE)
so_l <- lapply(so_l, RunUMAP, dims = 1:90, verbose = FALSE)
so_l <- lapply(so_l, RunTSNE, dims = 1:30, verbose = FALSE)

reduction <- "umap"

options(repr.plot.width = 25, repr.plot.height = 20)

plot <- function(so) {
    
    dplot_1 <- DimPlot(so, reduction = reduction, group.by = "RNA_snn_res.0.8", label = TRUE) & 
        ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
        theme(aspect.ratio = 1, legend.position = "none") & 
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
    
    dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) & 
        ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
        scale_color_manual(values = color$tissue, na.value = "dark gray") & 
        theme(aspect.ratio = 1, legend.position = "bottom") & 
        guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
    
    dplot_3 <- DimPlot(so, reduction = reduction, group.by = "sample_rep", label = FALSE) & 
        ggtitle("Replicate") & xlab("UMAP 1") & ylab("UMAP 2") &
        theme(aspect.ratio = 1, legend.position = "bottom") & 
        guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))

    dplot_4 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) &
        theme(aspect.ratio = 1, legend.position = "bottom") &
        scale_color_manual(values = color$main_labels, na.value = "dark gray") &
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

    dplot_5 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) &
        theme(aspect.ratio = 1, legend.position = "bottom") &
        scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
    
    fplot_1 <- FeaturePlot(so, reduction = reduction, features = "Cd34") & theme(aspect.ratio = 1)
    
    fplot_2 <- FeaturePlot(so, reduction = reduction, features = "Kit") & theme(aspect.ratio = 1)
    
    fplot_3 <- FeaturePlot(so, reduction = reduction, features = "Gata2") & theme(aspect.ratio = 1)
    
    fplot_4 <- FeaturePlot(so, reduction = reduction, features = "Gata1") & theme(aspect.ratio = 1)
    
    fplot_5 <- FeaturePlot(so, reduction = reduction, features = "Tfrc") & theme(aspect.ratio = 1) & ggtitle("Tfrc (CD71)")
    
    fplot_6 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
    
    fplot_7 <- FeaturePlot(so, reduction = reduction, features = "msTcr_RNA1") & theme(aspect.ratio = 1)
    
    fplot_8 <- FeaturePlot(so, reduction = reduction, features = "msIgkc_RNA1") & theme(aspect.ratio = 1)
    
    fplot_9 <- FeaturePlot(so, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
    
    fplot_10 <- FeaturePlot(so, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
    
    fplot_11 <- FeaturePlot(so, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)

    dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + fplot_7 + fplot_8 + fplot_9 + fplot_10 + fplot_11 + plot_layout(ncol = 4) 
    
}

lapply(so_l, plot)

$CpG

$NaCl

Remove cells by cluster expressing marker genes

so_cpg <- so_l[[1]]
so_nacl <- so_l[[2]]

B-cell cluster

CpG: Cluster 17
NaCl: Cluster 17

so_cpg <- subset(so_cpg, subset = RNA_snn_res.0.8 != 17)
so_nacl <- subset(so_nacl, subset = RNA_snn_res.0.8 != 17)

T-Cell cluster

CpG: None
NaCl: Cluster 18

so_nacl <- subset(so_nacl, subset = RNA_snn_res.0.8 != 18)

FACS sort overlap

options(repr.plot.width = 10, repr.plot.height = 5)

tissue_cpg_cluster <- ggplot(so_cpg@meta.data, aes(x = seurat_clusters, fill = tissue)) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$tissue) +
  ggtitle("CpG FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom")

tissue_nacl_cluster <- ggplot(so_nacl@meta.data, aes(x = seurat_clusters, fill = tissue)) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$tissue) +
  ggtitle("NaCl FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom")

tissue_cpg_cluster + tissue_nacl_cluster + plot_layout(ncol = 2)

Sort out cells which overlap FACS sort in low frequncy (N<=2)

clov: cluster overlap
clpu: cluster pure
clmi: cluster mixed

cpg_clov <- dplyr::group_by(so_cpg@meta.data, RNA_snn_res.0.8) %>% count(tissue)

# Get all pure cluster with no FACS mixture 
cpg_clpu <- cpg_clov[!cpg_clov$RNA_snn_res.0.8 %in% cpg_clov$RNA_snn_res.0.8[duplicated(cpg_clov$RNA_snn_res.0.8)], ] 

# Get all mixed cluster
cpg_clmi <- cpg_clov[!cpg_clov$RNA_snn_res.0.8 %in% cpg_clpu$RNA_snn_res.0.8, ] 

# Remove cells with low frequncy mixed cluster from Seurat object 
cpg_clmi_lf <- cpg_clmi[cpg_clmi$n <= 2, ] # Get all low frequency mixed cluster
so_cpg$clmi_lf <- ifelse(paste0(so_cpg$`RNA_snn_res.0.8`, so_cpg$tissue) %in% paste0(cpg_clmi_lf$RNA_snn_res.0.8, cpg_clmi_lf$tissue), TRUE, FALSE)
so_cpg <- subset(so_cpg, subset = clmi_lf == FALSE)

# Get all mixed cluster with high frequency for myeloid dominant and progenitor dominant 
cpg_clmi_hf <- cpg_clmi[!cpg_clmi$RNA_snn_res.0.8 %in%  cpg_clmi_lf$RNA_snn_res.0.8, ]
cpg_clmi_hf <- dplyr::group_by(cpg_clmi_hf, RNA_snn_res.0.8) %>% dplyr::filter(n == min(n))

cpg_clmi_hf_m <- cpg_clmi_hf[cpg_clmi_hf$tissue == "Myeloid", ]
cpg_clmi_hf_p <- cpg_clmi_hf[cpg_clmi_hf$tissue == "Progenitor", ]

nacl_clov <- dplyr::group_by(so_nacl@meta.data, RNA_snn_res.0.8) %>% count(tissue)

# Get all pure cluster with no FACS mixture 
nacl_clpu <- nacl_clov[!nacl_clov$RNA_snn_res.0.8 %in% nacl_clov$RNA_snn_res.0.8[duplicated(nacl_clov$RNA_snn_res.0.8)], ] 

# Get all mixed cluster
nacl_clmi <- nacl_clov[!nacl_clov$RNA_snn_res.0.8 %in% nacl_clpu$RNA_snn_res.0.8, ] 

# Remove cells with low frequncy mixed cluster from Seurat object 
nacl_clmi_lf <- nacl_clmi[nacl_clmi$n <= 2, ] # Get all low frequency mixed cluster
so_nacl$clmi_lf <- ifelse(paste0(so_nacl$`RNA_snn_res.0.8`, so_nacl$tissue) %in% paste0(nacl_clmi_lf$RNA_snn_res.0.8, nacl_clmi_lf$tissue), TRUE, FALSE)
so_nacl <- subset(so_nacl, subset = clmi_lf == FALSE)

# Get all mixed cluster with high frequency for myeloid dominant and progenitor dominant 
nacl_clmi_hf <- nacl_clmi[!nacl_clmi$RNA_snn_res.0.8 %in%  nacl_clmi_lf$RNA_snn_res.0.8, ]
nacl_clmi_hf <- dplyr::group_by(nacl_clmi_hf, RNA_snn_res.0.8) %>% dplyr::filter(n == min(n))

nacl_clmi_hf_m <- nacl_clmi_hf[nacl_clmi_hf$tissue == "Myeloid", ]
nacl_clmi_hf_p <- nacl_clmi_hf[nacl_clmi_hf$tissue == "Progenitor", ]

Find DEG in cluster with mixed Myeloid / Progenitor population

find_tissue_markers <- function(so, cluster, low_freq, high_freq, min_cells_group = 3) {
    
    so <- subset(so, subset = RNA_snn_res.0.8 %in% cluster)
    
    # Set idents to cluster tissue combination 
    so$cluster_tissue <- paste0(so$`RNA_snn_res.0.8`, "_", so$tissue)
    Idents(so) <- "cluster_tissue"
    
    # List to store results from looping through cluster
    deg_l <- list()
    
    for(cluster_idx in unique(so$`RNA_snn_res.0.8`)) {
        # Check if enough cells are present in both idents 
        so_cluster_idx <- subset(so, RNA_snn_res.0.8 == cluster_idx)
        sample_size_check <- all(table(so_cluster_idx$tissue) > min_cells_group) & length(unique(so_cluster_idx$tissue)) == 2
        # DEG
        if(sample_size_check) {
            deg <- FindMarkers(so, ident.1 = paste0(cluster_idx, "_", low_freq), ident.2 = paste0(cluster_idx, "_", high_freq), min.cells.group = min_cells_group, only.pos = TRUE, verbose = TRUE)
            # Filter deg list
            if(nrow(deg) > 0) {
                deg <- deg[grepl("^(mt-|Rp(s|l)|Gm\\d)", rownames(deg)) == FALSE, ]
                deg <- deg[!rownames(deg) %in% cc_kowalczyk$gene, ]  
            }
            deg_l[[cluster_idx]] <- deg

        } else {
            deg_l[[cluster_idx]] <- NULL
        }
    }
    
    return(deg_l)
}

deg_cpg_clmi_hf_m <- find_tissue_markers(so_cpg, cpg_clmi_hf_m$RNA_snn_res.0.8, low_freq = "Myeloid", high_freq = "Progenitor")
deg_cpg_clmi_hf_p <- find_tissue_markers(so_cpg, cpg_clmi_hf_p$RNA_snn_res.0.8, low_freq = "Progenitor", high_freq = "Myeloid")

deg_nacl_clmi_hf_m <- find_tissue_markers(so_nacl, nacl_clmi_hf_m$RNA_snn_res.0.8, low_freq = "Myeloid", high_freq = "Progenitor")
# deg_nacl_clmi_hf_p <- find_tissue_markers(so_nacl, nacl_clmi_hf_p$RNA_snn_res.0.8, low_freq = "Progenitor", high_freq = "Myeloid")

top_hm <- function(cluster_idx, so, deg, top = 25) {
    
    # Filter deg list
    deg <- deg[[cluster_idx]]
    
    # Get up and downreagulated genes
    deg_up <- deg[sign(deg$avg_log2FC) == 1, ]
    deg_down <- deg[sign(deg$avg_log2FC) == -1, ]

    # Select top hits by adjusted p value
    deg_up <- arrange(deg_up, p_val_adj)[1:top, ]
    deg_down <- arrange(deg_down, p_val_adj)[1:top, ]

    # Combine top hits 
    deg <- rbind(deg_up, deg_down)

    # Extract normalized count matrix for top hits
    so <- subset(so, subset = RNA_snn_res.0.8 == cluster_idx)
    cnt_norm <- GetAssayData(so, assay = "RNA", slot = "data")
    cnt_norm <- cnt_norm[rownames(cnt_norm) %in% rownames(deg), ]

    # Order count matrix by tissue 
    meta <- so@meta.data
    meta <- arrange(meta, tissue)
    cnt_norm <- cnt_norm[, rownames(meta)]
    
    # Compute mean per group
    colnames(cnt_norm) <- meta$tissue
    cnt_norm <- data.frame(
        Myeloid = rowMeans(cnt_norm[, grepl("Myeloid", colnames(cnt_norm), fixed = TRUE)]), 
        Progenitor = rowMeans(cnt_norm[, grepl("Progenitor", colnames(cnt_norm), fixed = TRUE)]))

    # ComplexHeatmap color ramp
    color_range <- max(abs(cnt_norm))
    color_break <- seq(0, color_range, 0.01)
    color_ramp <- viridis(length(color_break), option = "magma")
    
    # Top annotaion 
    annoation_col <- data.frame("FACS" = colnames(cnt_norm))
    rownames(annoation_col) <- colnames(cnt_norm)
    annotation_colors <- list("FACS" = color$tissue)
    
    # Heat map 
    hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
        mat           = as.matrix(cnt_norm),
        main          = paste("Cluster", cluster_idx),
        fontsize_row  =  10,
        scale         = "none",
        cluster_rows  = TRUE,
        cluster_cols  = FALSE,
        cellwidth      = 12, 
        cellheight     = 12, 
        clustering_distance_rows = "euclidean",
        clustering_distance_cols = "euclidean",
        clustering_method        = "complete",
        show_row_dend = FALSE, 
        show_rownames = TRUE,
        show_colnames = FALSE,
        annotation_col = annoation_col, 
        annotation_colors = annotation_colors, 
        color         = color_ramp,
        breaks        = color_break, 
        border_color  = NA, 
        use_raster = TRUE)))
    
    return(hm)
}

# CpG
hm_cpg_clmi_hf_m <- lapply(names(deg_cpg_clmi_hf_m), top_hm, so = so_cpg, deg = deg_cpg_clmi_hf_m)
hm_cpg_clmi_hf_p <- lapply(names(deg_cpg_clmi_hf_p), top_hm, so = so_cpg, deg = deg_cpg_clmi_hf_p)
# NaCl
hm_nacl_clmi_hf_m <- lapply(names(deg_nacl_clmi_hf_m), top_hm, so = so_nacl, deg = deg_nacl_clmi_hf_m)

Cpg Myeloid cells mixed in Progenitor dominant cluster

options(repr.plot.width = 20, repr.plot.height = 10)
gridExtra::arrangeGrob(grobs = hm_cpg_clmi_hf_m, nrow = 2) %>% grid::grid.draw()

so_cpg$clmi_hf <- ifelse(paste0(so_cpg$`RNA_snn_res.0.8`, so_cpg$tissue) %in% paste0(c(2, 8, 1, 4, 0, 6, 3, 12, 7), "Myeloid"), TRUE, FALSE)
so_cpg <- subset(so_cpg, subset = clmi_hf == FALSE)

Cpg Progenitor cells mixed in Myeloid dominant cluster

Keep all. It looks like a mixture of myeloid and MEP progenitors

options(repr.plot.width = 20, repr.plot.height = 5)
gridExtra::arrangeGrob(grobs = hm_cpg_clmi_hf_p, nrow = 1) %>% grid::grid.draw()

NaCl Myeloid cells mixed in Progenitor dominant cluster

Only keep cluster 14 because it looks like a true myeloid progenitor cluster

options(repr.plot.width = 20, repr.plot.height = 10)
gridExtra::arrangeGrob(grobs = hm_nacl_clmi_hf_m, nrow = 2) %>% grid::grid.draw()

so_nacl$clmi_hf <- ifelse(paste0(so_nacl$`RNA_snn_res.0.8`, so_nacl$tissue) %in% paste0(c(5, 7, 12, 6, 10, 0), "Myeloid"), TRUE, FALSE)
so_nacl <- subset(so_nacl, subset = clmi_hf == FALSE)

Recompute dimensional reduction

so_l <- list(so_cpg, so_nacl)

so_l <- lapply(so_l, FindVariableFeatures, nfeatures = 2000)

so_l <- lapply(so_l, variable_features_reduce, genes = cc_kowalczyk$gene)

so_l <- lapply(so_l, ScaleData, vars.to.regress = c("nCount_RNA"))

Regressing out nCount_RNA

Centering and scaling data matrix

Regressing out nCount_RNA

Centering and scaling data matrix

options(warn = -1)
so_l <- lapply(so_l, RunPCA, npcs = 100, verbose = FALSE)
so_l <- lapply(so_l, FindNeighbors, dims = 1:30, verbose = FALSE)
so_l <- lapply(so_l, FindClusters, verbose = FALSE)
so_l <- lapply(so_l, RunUMAP, dims = 1:90, verbose = FALSE)
so_l <- lapply(so_l, RunTSNE, dims = 1:30, verbose = FALSE)

reduction <- "umap"

options(repr.plot.width = 25, repr.plot.height = 20)

plot <- function(so) {
    
    dplot_1 <- DimPlot(so, reduction = reduction, group.by = "RNA_snn_res.0.8", label = TRUE) & 
        ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
        theme(aspect.ratio = 1, legend.position = "none") & 
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
    
    dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) & 
        ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
        scale_color_manual(values = color$tissue, na.value = "dark gray") & 
        theme(aspect.ratio = 1, legend.position = "bottom") & 
        guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))
    
    dplot_3 <- DimPlot(so, reduction = reduction, group.by = "sample_rep", label = FALSE) & 
        ggtitle("Replicate") & xlab("UMAP 1") & ylab("UMAP 2") &
        theme(aspect.ratio = 1, legend.position = "bottom") & 
        guides(color = guide_legend(ncol = 1, override.aes = list(size = 2)))

    dplot_4 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) &
        theme(aspect.ratio = 1, legend.position = "bottom") &
        scale_color_manual(values = color$main_labels, na.value = "dark gray") &
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

    dplot_5 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) &
        theme(aspect.ratio = 1, legend.position = "bottom") &
        scale_color_manual(values = color$fine_labels, na.value = "dark gray") &
        guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
    
    fplot_1 <- FeaturePlot(so, reduction = reduction, features = "Cd34") & theme(aspect.ratio = 1)
    
    fplot_2 <- FeaturePlot(so, reduction = reduction, features = "Kit") & theme(aspect.ratio = 1)
    
    fplot_3 <- FeaturePlot(so, reduction = reduction, features = "Gata2") & theme(aspect.ratio = 1)
    
    fplot_4 <- FeaturePlot(so, reduction = reduction, features = "Gata1") & theme(aspect.ratio = 1)
    
    fplot_5 <- FeaturePlot(so, reduction = reduction, features = "Tfrc") & theme(aspect.ratio = 1) & ggtitle("Tfrc (CD71)")
    
    fplot_6 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)
    
    fplot_7 <- FeaturePlot(so, reduction = reduction, features = "msTcr_RNA1") & theme(aspect.ratio = 1)
    
    fplot_8 <- FeaturePlot(so, reduction = reduction, features = "msIgkc_RNA1") & theme(aspect.ratio = 1)
    
    fplot_9 <- FeaturePlot(so, reduction = reduction, features = "nCount_RNA") & theme(aspect.ratio = 1)
    
    fplot_10 <- FeaturePlot(so, reduction = reduction, features = "nFeature_RNA") & theme(aspect.ratio = 1)
    
    fplot_11 <- FeaturePlot(so, reduction = reduction, features = "pMt_RNA") & theme(aspect.ratio = 1)

    dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + fplot_6 + fplot_7 + fplot_8 + fplot_9 + fplot_10 + fplot_11 + plot_layout(ncol = 4) 
    
}

lapply(so_l, plot)

[[1]]

[[2]]

options(repr.plot.width = 10, repr.plot.height = 5)

tissue_cpg_cluster <- ggplot(so_l[[1]]@meta.data, aes(x = seurat_clusters, fill = tissue)) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$tissue) +
  ggtitle("CpG FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom")

tissue_nacl_cluster <- ggplot(so_l[[2]]@meta.data, aes(x = seurat_clusters, fill = tissue)) + 
  geom_bar(stat = "count", position = "fill") + 
  scale_fill_manual(values = color$tissue) +
  ggtitle("NaCl FACS sort frequency") + xlab("Cluster") + ylab("Cell frequency") + 
  theme(legend.position = "bottom")

tissue_cpg_cluster + tissue_nacl_cluster + plot_layout(ncol = 2)

Save output

so <- merge(so_l[[1]], so_l[[2]])

saveRDS(so, so_out_file)

TEST

so <- NormalizeData(
  object               = so, 
  assay                = "RNA", 
  normalization.method = "LogNormalize", 
  scale.factor         = 10000
  )

so <- FindVariableFeatures(
  object           = so, 
  assay            = "RNA", 
  selection.method = "vst", 
  nfeatures        = 3000
  )

variable_features <- VariableFeatures(so)
VariableFeatures(so) <- variable_features[!variable_features %in% c(cc_kowalczyk$gene)]

so <- ScaleData(
  object          = so, 
  assay           = "RNA", 
  vars.to.regress = c("nCount_RNA"), 
  do.scale        = TRUE, 
  do.center       = TRUE
  )

Regressing out nCount_RNA

Centering and scaling data matrix

source("bin/seurat_function.R")
so <- dim_clust(
    # Seurat 
    so = so, 
    assay = "RNA",
    # PCA features
    blacklist_genes = NULL, 
    # Dim reduction
    dims_pca  = 100,  # Default 20 - RunPCA dims
    dims_umap = 90,  # Default NULL - RunUMAP dims 
    dims_tsne = 30,  # Default 1:5
    min_dist  = 0.3,  # Default 0.3  - RunUMAP dmin.ist - controls how tightly the embedding
    # Cluster 
    dims_cluster = 10,  # Default 1:10 - FindNeighbors dims  
    cluster_res  = 0.8 # Default 0.8 - FindClusters resoluton - above(below) 1.0 to obtain larger(smaller) number of communities)
)

PC_ 1 
Positive:  Car2, Prdx2, Blvrb, Glrx5, Rps2, Hbb-bs, Hba-a1, Cd24a, Hbb-bt, Cpox 
	   Ctse, Rhd, Car1, Mgst3, Hba-a2, Cldn13, Odc1, Ybx3, Ermap, Stard10 
	   Hemgn, Gypa, C1qtnf12, Klf1, Rpl14, Hdgf, Hsp90aa1, Aqp1, Ubac1, Smim1 
Negative:  Tyrobp, Tmsb4x, Ptprc, Laptm5, Cd52, Fcer1g, Cyba, Ctss, Coro1a, Cd74 
	   Lcp1, Lsp1, Cst3, Sh3bgrl3, Napsa, Arpc1b, Gm2a, Spi1, Lst1, Gngt2 
	   Rgs2, Sat1, Psap, Fyb, Cybb, H2-Aa, Pou2f2, Pld4, Actb, H2-Ab1 
PC_ 2 
Positive:  Fth1, Ftl1, Hba-a2, Hbb-bt, Hba-a1, Hbb-bs, Ctsb, Gypa, Hebp1, Slc4a1 
	   Mgst3, Ctse, Hemgn, Blvrb, Cd24a, Rhd, Slc25a37, Cldn13, Hagh, Acp5 
	   Cpox, Tfrc, Alas2, Ubac1, Prdx2, Isg20, C1qc, Ctsd, Vcam1, C1qa 
Negative:  Myb, S100a10, Npm1, Tagln2, Gm15915, Ptprcap, Ldha, Cmtm7, Mif, Srm 
	   H2afy, Sox4, Phgdh, Nop10, Ybx1, Npm3, Actg1, Tmsb10, Ppp1r14b, Hsp90ab1 
	   Mt1, Prdx6, Fbl, Hspd1, Hspa8, Srgn, Sptssa, Tomm5, Slc25a4, Psmb8 
PC_ 3 
Positive:  Vcam1, C1qa, C1qc, C1qb, Fcna, Cd5l, Cadm1, Mertk, Igf1, Mrc1 
	   Axl, Frmd4b, Itgad, Maf, Actn1, Id2, Slc40a1, Spic, Ccr3, Dnase1l3 
	   Mafb, Hpgds, Asb2, Ptgs1, Cfp, Kcna2, Epb41l3, Slpi, Hmox1, Aif1 
Negative:  Ace, Cebpb, Itgal, Msrb1, Lyz2, Cyp4f18, Cx3cr1, Spn, Ceacam1, Cd300ld 
	   Tppp3, Smpdl3b, Nupr1, Hp, Klf4, Pglyrp1, Smpdl3a, Nr4a1, Eno3, Rbpms 
	   Fabp4, Ear2, Stk10, Ifitm2, Klf2, Tnfrsf1b, S1pr5, Il17ra, Rras, Dusp16 
PC_ 4 
Positive:  Vcam1, C1qc, C1qa, Fcna, Mertk, Cd5l, C1qb, Igf1, Axl, Mafb 
	   Mrc1, Maf, Itgad, Slc40a1, Hmox1, Ccr3, Apoe, Slpi, Spic, Ptgs1 
	   Epb41l3, Trf, Selenop, Kcna2, Kcnj10, Pla2g15, Clec4n, Sema6d, Timp2, Gfra2 
Negative:  Ciita, Ccnd1, H2-Eb1, H2-Oa, Dnase1l3, H2-Aa, H2-Ab1, Slamf7, H2-DMb2, Ppp1r14a 
	   Traf1, Lsp1, Tbc1d4, S100a11, Ffar2, Siglecg, Wdfy4, Jaml, Adam23, Cd74 
	   Avpi1, Rogdi, Lyst, Flt3, Ccl5, S100a4, H2-Ob, Cd52, Gsn, Cd7 
PC_ 5 
Positive:  Xcr1, Naaa, Rab7b, Irf8, Cd8a, Clec9a, Ppt1, P2ry14, Tlr3, Cadm1 
	   A530099J19Rik, AC163354.1, Ifi205, Mpeg1, Pianp, Pdia5, Itgae, Gcsam, Olfm1, Crip1 
	   Sept3, Lgals3, Tlr11, Aif1, Dbn1, Btla, Cxcl9, Slc8b1, Ece1, Cst3 
Negative:  Ccl5, S100a4, Ccnd1, Relb, Ffar2, Siglecg, Dtx1, Ltb, H2-DMb2, Il4i1 
	   Plxnc1, Tbc1d4, Cd7, Cyp4f16, Tmem176b, Traf1, Ppp1r14a, Tmem176a, H2-Eb2, Ryr3 
	   Man1a, Adam23, Lyst, Bcl2a1b, Avpi1, Icosl, Tctex1d2, Arap2, St8sia6, Ube2e2 

Computing nearest neighbor graph

Computing SNN

Computing nearest neighbors

Only one graph name supplied, storing nearest-neighbor graph only

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 21646
Number of edges: 732713

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9031
Number of communities: 22
Elapsed time: 3 seconds

17:33:05 UMAP embedding parameters a = 0.9922 b = 1.112

17:33:06 Commencing smooth kNN distance calibration using 1 thread

17:33:07 Initializing from normalized Laplacian + noise

17:33:08 Commencing optimization for 200 epochs, with 567640 positive edges

17:33:31 Optimization finished

17:33:31 UMAP embedding parameters a = 0.9922 b = 1.112

17:33:31 Read 21646 rows and found 90 numeric columns

17:33:31 Using Annoy for neighbor search, n_neighbors = 20

17:33:31 Building Annoy index with metric = euclidean, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

17:33:33 Writing NN index file to temp file /tmp/RtmpMWeusO/file31487b7b8bd4

17:33:33 Searching Annoy index using 1 thread, search_k = 2000

17:33:37 Annoy recall = 100%

17:33:38 Commencing smooth kNN distance calibration using 1 thread

17:33:39 Initializing from normalized Laplacian + noise

17:33:40 Commencing optimization for 200 epochs, with 682638 positive edges

17:34:05 Optimization finished

reduction <- "rna_umap_nno"

dplot_1 <- DimPlot(so, reduction = reduction, group.by = "rna_snn_res.0.8", label = TRUE) & 
    theme(aspect.ratio = 1, legend.position = "none")
dplot_2 <- DimPlot(so, reduction = reduction, group.by = "main_labels", label = FALSE) & 
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$main_labels, na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_3 <- DimPlot(so, reduction = reduction, group.by = "fine_labels", label = FALSE) & 
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$fine_labels, na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_4 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) & 
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$tissue, na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(so, reduction = reduction, group.by = "treatment", label = FALSE) & 
    theme(aspect.ratio = 1, legend.position = "bottom") & 
    scale_color_manual(values = color$treatment, na.value = "dark gray") & 
    guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))

fplot_1 <- FeaturePlot(so, reduction = reduction, features = "pHb_RNA") & theme(aspect.ratio = 1)

options(repr.plot.width = 15, repr.plot.height = 10)
dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(ncol = 3)
dplot

Session info

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ComplexHeatmap_2.8.0 viridis_0.6.1        viridisLite_0.4.0   
 [4] gridExtra_2.3        patchwork_1.1.1      RColorBrewer_1.1-2  
 [7] ggplot2_3.3.3        dplyr_1.0.6          SeuratObject_4.0.1  
[10] Seurat_4.0.2        

loaded via a namespace (and not attached):
  [1] Rtsne_0.15            colorspace_2.0-1      rjson_0.2.20         
  [4] deldir_0.2-10         ellipsis_0.3.2        ggridges_0.5.3       
  [7] circlize_0.4.13       IRdisplay_1.0         GlobalOptions_0.1.2  
 [10] base64enc_0.1-3       clue_0.3-59           spatstat.data_2.1-0  
 [13] farver_2.1.0          leiden_0.3.8          listenv_0.8.0        
 [16] ggrepel_0.9.1         RSpectra_0.16-0       fansi_0.5.0          
 [19] codetools_0.2-18      splines_4.1.0         doParallel_1.0.16    
 [22] polyclip_1.10-0       IRkernel_1.2          jsonlite_1.7.2       
 [25] Cairo_1.5-12.2        ica_1.0-2             cluster_2.1.2        
 [28] png_0.1-7             uwot_0.1.10           shiny_1.6.0          
 [31] sctransform_0.3.2     spatstat.sparse_2.0-0 compiler_4.1.0       
 [34] httr_1.4.2            assertthat_0.2.1      Matrix_1.3-4         
 [37] fastmap_1.1.0         lazyeval_0.2.2        limma_3.48.1         
 [40] later_1.2.0           htmltools_0.5.1.1     tools_4.1.0          
 [43] igraph_1.2.6          gtable_0.3.0          glue_1.4.2           
 [46] RANN_2.6.1            reshape2_1.4.4        Rcpp_1.0.6           
 [49] scattermore_0.7       vctrs_0.3.8           nlme_3.1-152         
 [52] iterators_1.0.13      lmtest_0.9-38         stringr_1.4.0        
 [55] globals_0.14.0        mime_0.10             miniUI_0.1.1.1       
 [58] lifecycle_1.0.0       irlba_2.3.3           goftest_1.2-2        
 [61] future_1.21.0         MASS_7.3-54           zoo_1.8-9            
 [64] scales_1.1.1          spatstat.core_2.1-2   promises_1.2.0.1     
 [67] spatstat.utils_2.1-0  parallel_4.1.0        reticulate_1.20      
 [70] pbapply_1.4-3         rpart_4.1-15          stringi_1.6.2        
 [73] S4Vectors_0.30.0      foreach_1.5.1         BiocGenerics_0.38.0  
 [76] shape_1.4.6           repr_1.1.3            rlang_0.4.11         
 [79] pkgconfig_2.0.3       matrixStats_0.59.0    evaluate_0.14        
 [82] lattice_0.20-44       ROCR_1.0-11           purrr_0.3.4          
 [85] tensor_1.5            labeling_0.4.2        htmlwidgets_1.5.3    
 [88] cowplot_1.1.1         tidyselect_1.1.1      parallelly_1.25.0    
 [91] RcppAnnoy_0.0.18      plyr_1.8.6            magrittr_2.0.1       
 [94] R6_2.5.0              magick_2.7.2          IRanges_2.26.0       
 [97] generics_0.1.0        pbdZMQ_0.3-5          DBI_1.1.1            
[100] pillar_1.6.1          withr_2.4.2           mgcv_1.8-36          
[103] fitdistrplus_1.1-5    survival_3.2-11       abind_1.4-5          
[106] tibble_3.1.2          future.apply_1.7.0    crayon_1.4.1         
[109] uuid_0.1-4            KernSmooth_2.23-20    utf8_1.2.1           
[112] spatstat.geom_2.1-0   plotly_4.9.3          GetoptLong_1.0.5     
[115] data.table_1.14.0     digest_0.6.27         xtable_1.8-4         
[118] tidyr_1.1.3           httpuv_1.6.1          stats4_4.1.0         
[121] munsell_0.5.0